/*
 * ReactionSystem.cpp
 *
 *  Created on: 23 Jul 2011
 *      Author: Allan
 */

#include "ReactionSystem.h"

// C++ includes
#include <iomanip>
#include <set>

// Boost includes
#include <boost/foreach.hpp>

// Eigen includes
#include <Eigen/Core>
using namespace Eigen;

// GeoReact includes
#include "Utils.h"

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/// Basis Swap Implementation todo finish this implementation
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//void ReactionSystem::RestructureReactions(const vector<string>& newPrimarySpecies, const vector<string>& newSecondarySpecies)
//{
//	unsigned Nj = products.size();
//	unsigned Ni = reactants.size();
//	
//	MatrixXd M = MatrixXd::Zero(Nj+Ni, Nj);
//	
//	M.block( 0, 0, Nj, Nj) = MatrixXd::Identity(Nj, Nj);
//	M.block(Nj, 0, Ni, Nj) = vij;
//	
//	// The primary-swap matrix (theta)
//	MatrixXd theta = MatrixXd::Zero(Nj, Nj);
//	
//	// Loop over all the species in (newPrimarySpecies) and set the primary-swap matrix (theta) 
//	for(unsigned row = 0; row < Nj; ++row)
//	{
//		// Determine the local indexes of the current primary species in the set of the actual primary and secondary species, (primarySpecies) and (secondarySpecies) respectively (note that only one of the indexes below will be within range) 
//		const unsigned j = std::find(products.begin(), products.end(), newPrimarySpecies[row]) - products.begin();
//		const unsigned i = std::find(reactants.begin(), reactants.end(), newPrimarySpecies[row]) - reactants.begin();
//		
//		// Check if the current new primary species (primary) is in the set of primary species (primarySpecies) and set the row of (theta) accordingly 
//		if(j < Nj) theta(row, j) = 1.0;
//		
//		// Check if the current new primary species (primary) is in the set of secondary species (secondarySpecies) and set the row of (theta) accordingly 
//		if(i < Ni) theta.row(row) = vij.row(i);
//	}
//	
//	// Transform the matrix M to reflect the changes in the set of primary species (M' = M * inv(theta))
//	MatrixXd thetaInv = FullPivLU<MatrixXd>(theta).inverse();
//	
//	M *= thetaInv;
//	
//	// Update the stoichiometric matrix vij
//	for(unsigned row = 0; row < Ni; ++row)
//	{
//		// Determine the local indexes of the current primary species in the set of the actual primary and secondary species, (primarySpecies) and (secondarySpecies) respectively (note that only one of the indexes below will be within range) 
//		const unsigned j = std::find(products.begin(), products.end(), newSecondarySpecies[row]) - products.begin();
//		const unsigned i = std::find(reactants.begin(), reactants.end(), newSecondarySpecies[row]) - reactants.begin();
//		
//		vij.row(row) = (j < Nj) ? M.row(j) : M.row(i + Nj);
//	}
//	
//	// Create the new reactions
//	vector<Reaction> newReactions(Ni);
//	
//	for(unsigned i = 0; i < Ni; ++i)
//	{
//		vector<string> products;
//		vector<double> stoichiometries;
//		
//		const double tol = 1.0E-6;
//		
//		for(unsigned j = 0; j < Nj; ++j)
//		{
//			if(std::abs(vij(i, j)) < tol)
//			{
//				products.push_back(newPrimarySpecies[j]);
//				stoichiometries.push_back(vij(i, j));
//			}
//		}
//		
//		newReactions[i].SetEquation(newSecondarySpecies[i], products, stoichiometries);
//	}
//	
//	
//}

ReactionSystem::ReactionSystem()
{}

ReactionSystem::ReactionSystem(const vector<Reaction>& reactions) :
reactions(reactions)
{}

const vector<Reaction>& ReactionSystem::GetReactions() const
{
	return reactions;
}

const vector<string> ReactionSystem::Reactants() const
{
	// The reactants of the reactions
	vector<string> reactants;
	
	// Loop over all the reactions
	BOOST_FOREACH(const Reaction& reaction, reactions)
		reactants.push_back(reaction.GetReactant().name);
	
	return reactants;
}

const vector<string> ReactionSystem::Products() const
{
	// The unique products of the reactions
	set<string> products;
	
	// Loop over all the reactions and insert their products into the set of products
	BOOST_FOREACH(const Reaction& reaction, reactions)
		BOOST_FOREACH(const ReactionSpecies& product, reaction.GetProducts())
			products.insert(product.name);

	return vector<string>(products.begin(), products.end());
}

const double ReactionSystem::Stoichiometry(const string& reactant, const string& product) const
{
	// An auxiliar function to determine if the reactant of a given reaction is (reactant)
	auto checkReactant = [&](const Reaction& reaction) { return reaction.GetReactant().name == reactant; };
	
	// The reaction whose reactant is (reactant)
	const Reaction& reaction = *find_if(reactions.begin(), reactions.end(), checkReactant);
	
	// The products of the reaction
	const vector<ReactionSpecies>& products = reaction.GetProducts();
	
	// Loop over all the products in (products) and return the stoichiometry of the product that matches the name of the given one (product)
	BOOST_FOREACH(const ReactionSpecies& p, products) 
		if(p.name == product) return p.stoichiometry;
	
	// Return zero, since (product) does not compose the given (reactant)
	return 0.0;
}

const MatrixXd ReactionSystem::StoichiometricMatrix(const vector<string>& reactants, const vector<string>& products) const
{
	// The number of reactants (Nr) and products (Np)
	const unsigned Nr = reactants.size();
	const unsigned Np = products.size();

	// The stoichiometric matrix for the list of (reactants) and (products)
	MatrixXd v = MatrixXd::Zero(Nr, Np);
	
	// Create the stoichiometric matrix
	for(unsigned r = 0; r < Nr; ++r) 
		for(unsigned p = 0; p < Np; ++p)
			v(r, p) = Stoichiometry(reactants[r], products[p]);
	
	return v;
}

const VectorXd ReactionSystem::EquilibriumConstants(double T, double P) const
{
	VectorXd Ki(reactions.size());
	
	for(int i = 0; i < Ki.rows(); ++i)
		Ki[i] = reactions[i].GetEquilibriumConstant()(T, P);
	
	return Ki;
}

const VectorXd ReactionSystem::ReactionQuotients(const VectorXd& a) const
{
	VectorXd Qi(reactions.size());
	
	for(int i = 0; i < Qi.rows(); ++i)
		Qi[i] = reactions[i].ReactionQuotient(a);
	
	return Qi;
}

const VectorXd ReactionSystem::MassRates(double T, double P, const VectorXd& n, const VectorXd& a) const
{
	VectorXd Ri(reactions.size());
	
	for(int i = 0; i < Ri.rows(); ++i)
		Ri[i] = reactions[i].MassRate(T, P, n, a);
	
	return Ri;
}

void ReactionSystem::PrintState(double T, double P, const VectorXd& n, const VectorXd& a) const
{
	cout << setw(25) << left << "Reactants";
	cout << setw(25) << left << "log(Qr)";
	cout << setw(25) << left << "log(Kr)";
	cout << setw(25) << left << "log(Kr/Qr)" << endl << endl;
	
	BOOST_FOREACH(auto reaction, reactions)
	{
		const double Qr = reaction.ReactionQuotient(a);
		const double Kr = reaction.GetEquilibriumConstant()(T, P);
		
		cout << setw(25) << left << reaction.GetReactant().name;
		cout << setw(25) << left << std::log10(Qr);
		cout << setw(25) << left << std::log10(Kr);
		cout << setw(25) << left << std::log10(Kr/Qr) << endl;
	}
}

ostream& operator<<(ostream& out, const ReactionSystem& reactionSystem)
{
	out << "Reactions: " << endl;
	
	BOOST_FOREACH(auto reaction, reactionSystem.GetReactions())
		out << '\t' << reaction << endl;
	
	out << endl;
	
	out << "Reactants: " << endl << '\t';
		out << reactionSystem.Reactants() << endl << endl;
		
	out << "Products: " << endl << '\t';
		out << reactionSystem.Products() << endl << endl;
		
	out << "Stoichiometric Matrix: " << endl;
		out << reactionSystem.StoichiometricMatrix(reactionSystem.Reactants(), reactionSystem.Products()) << endl;
		
	return out;
}

